In [1]:
from IPython.core.display import Image
%matplotlib inline
print "Methanol molecule"
Image(filename='methanol.png')
Out[1]:
In [2]:
print "A galaxy at a lookback time of 7 billion years"
Image(filename='pks1830.jpg')
Out[2]:
In [3]:
from math import *
import numpy as np
from pylab import *
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.optimize import curve_fit
from linfit import linfit
from matplotlib import animation
from JSAnimation.IPython_display import display_animation
plt.rcParams['figure.figsize'] = 10, 6
# Linelist Nr 1:
name_list = "Linelist Nr 1"
transitions = [[22.2, 3],[142.2, -1.9],[106.4, 0],[46.0, 1.]] # [Tau [10^-3], Kmu]
# Some constants
sensitivity = 4.0*10**(-3) # [Jy]
i_bg = 2.0 # source intensity [Jy]
f_c = 0.38 # source covering factor
dmu = 5.0*10**(-6) # input dmu
fwhm_line = 20. # FWHM width [km/s]
sigma_line = fwhm_line/2.35 # sigma width [km/s]
res_bin = 1.8 # resolution bin [km/s]
sol = 299792.458 # c [km/s]
# Global functions
# Returns line strength in Jy, its expected SNR, and Kmu
def absline(line_list, i):
tau_line = line_list[i][0]*10**(-3)
i_line = (1-exp(-tau_line))*f_c*i_bg # line strength in [Jy]
snr = round(i_line/sensitivity,1) # snr = line_strength/sensitivity
Kmu = line_list[i][1]
return i_line, snr, Kmu
# Returns a gaussian profile
def gaus(x,a,x0,sigma):
return -a*np.exp(-(x-x0)**2/(2*sigma**2))
In [4]:
Nframes = 50
def animate(mu):
mu = (mu-Nframes/2)*10**-6
begin, end = -260,260
N = int(round((end - begin)/res_bin,0))
x = np.linspace(begin,end,N)
y=np.zeros_like(x)
for k in np.arange(len(transitions)):
i_line, snr, Kmu = absline(transitions, k)
center = -mu*Kmu*sol -195 + k*130
y += gaus(x,i_line,center,sigma_line)
l.set_data(x, y)
return l, dmu
fig, ax = plt.subplots()
ax = plt.axes(xlim=(-260,260), ylim=(-0.12, 0.02))
ax.tick_params(axis='both', which='major', labelsize=14)
ax.set_ylabel("Line strength", fontsize = 14)
ax.set_xlabel("Velocity [km/s]", fontsize = 14)
ax.set_title(r"Line positions depend uniquely on $\mu$", fontsize = 16, color="blue")
fig.tight_layout()
l, = ax.plot([],[])
begin, end = -260,260
N = int(round((end - begin)/res_bin,0))
x = np.linspace(begin,end,N)
y=np.zeros_like(x)
for k in np.arange(len(transitions)):
i_line, snr, Kmu = absline(transitions, k)
center = 0.0*Kmu*sol -195 + k*130
y += gaus(x,i_line,center,sigma_line)
ax.plot(x,y,'-k', label = "Measured in the lab")
lg = legend(loc = "lower right", prop={'size':16})
lg.draw_frame(False)
#print(animate(Nframes)[1])
anim = animation.FuncAnimation(fig, animate, frames=Nframes, interval=50, blit=True)
display_animation(anim, default_mode='once')
Out[4]:
In [171]:
# Number of spectra to be generated
nr_of_spectra =1000
output_dmu = []
output_dmu_err = []
for i in range(int(nr_of_spectra)):
# For each transition simulate a gaussian profile, add noise and fit
fitted = []
for k in np.arange(len(transitions)):
i_line, snr, Kmu = absline(transitions, k)
begin, end = -70, 70
N = round((end - begin)/res_bin,0)
x = np.linspace(begin,end,N)
center = -dmu*Kmu*sol
y = gaus(x,i_line,center,sigma_line) # returns a gaussian
yn = y + np.random.normal(0, sensitivity, size=len(x)) # noise added to gaussian
popt, pcov = curve_fit(gaus, x, yn) # fitting; popt = OPTimal Parameters for fit; COVariance matrix
#sigma = np.sqrt([pcov[0,0], pcov[1,1], pcov[2,2]]) # uncertainties from gaussian fitting
uncertainty_v = np.sqrt(pcov[1,1])
# Plot one realization to check if it looks ok
if i == 0:
#print i_line, sqrt(pcov[1,1]), uncertainty_v
fig, ax = plt.subplots()
ax.plot(x, yn, label=r"K$_{\mu}$ ="+str(Kmu)+" SNR="+str(snr))
ax.plot(x, gaus(x,popt[0],popt[1],popt[2]))
ax.tick_params(axis='both', which='major', labelsize=14)
legend = ax.legend(loc='lower right', shadow=True)
plt.title("Single realization", fontsize=14)
plt.xlabel("Velocity [km/s]", fontsize=14)
plt.ylabel("Line strength [Jy]", fontsize=14)
ax.set_ylim([-0.11,0.02])
#print popt, sigma # fitted parameters (intensity, centroid position, width) and their uncertainties
try:
fitted.append([Kmu, round(popt[1],4), round(uncertainty_v,4)])
except TypeError:
print "Type Error"
continue
# Form 3 arrays: Kmu, V/c, and (V/c)_error
x,y,sigmay = [], [], []
for j in np.arange(len(fitted)):
x.append(fitted[j][0])
y.append(fitted[j][1]/sol)
sigmay.append(fitted[j][2]/sol)
# Find dmu/mu
# Fit a linear model to the simulated data set; uses linfit with weighting
fit, cvm, info = linfit(x, y, sigmay, relsigma=False, return_all=True)
dfit = [np.sqrt(cvm[i,i]) for i in range(2)]
output_dmu_err.append([dfit[0]*np.sqrt(info.rchisq)])
output_dmu.append(-fit[0])
# Print some extra stuff
#print u"redchisq =", info.rchisq
#print u"slope =", fit[0], dfit[0]
#print "weighted slope uncertainty =", dfit[0]*np.sqrt(info.rchisq)
#print u"y-intercept =", fit[1], dfit[1]
In [173]:
# Print results
dmu_graph = plt.hist(output_dmu, bins = 20)
plt.legend(loc='upper left', shadow=True)
plt.ticklabel_format(style='sci', axis='x', scilimits=(0,0))
plt.title(r"$\Delta\mu/\mu$ values from %d artificial spectra" % (nr_of_spectra,), fontsize=16)
plt.xlabel(r"$\Delta\mu/\mu$", fontsize=18)
plt.ylabel("Frequency", fontsize=16)
plt.tick_params(axis='both', which='major', labelsize=14)
plt.savefig("linelist1.pdf", format='PDF')
print "Sensitivity:", str('%.1e' % sensitivity), "[Jy]"
print "Resolution bin:", str(res_bin), "[km/s]"
print "Source intensity:", str(i_bg), "[Jy]"
print "Covering factor:", str(f_c)
print "FWHM:", str(fwhm_line), "[km/s]"
print name_list
print "Number of methanol transitions:", len(transitions)
print "K_mu:", [transitions[i][1] for i in range(len(transitions))]
print "tau:", [transitions[i][0] for i in range(len(transitions))], "x10^-3"
print "---------------------------------"
(mu, sigma) = norm.fit(output_dmu)
plt.errorbar(x=mu, xerr=sigma, y=np.max(dmu_graph[0])/1.6, fmt='ro',markersize=12, linewidth=2, capsize=12, markeredgewidth=2)
print "Input dmu/mu =", '%.2e' % dmu
print "Output dmu/mu =", '%.2e' % mu, "+/-", '%.2e' %sigma
In [ ]: